﻿using System;
using System.Diagnostics;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using Nature.Numerics;
using Nature.Numerics.Basic;

namespace Nature
{
    public enum TemperatureChangeType
    {
        AtConstantPressure,
        AtConstantVolume,
        IsentropicCompression,
        ShockCompression
    }

    public enum PressureChangeType
    {
        AtConstantTemperature,
        AtConstantVolume,
        IsentropicCompression,
        ShockCompression
    }

    public enum DensityChangeType
    {
        AtConstantTemperature,
        AtConstantPressure,
        IsentropicCompression,
        ShockCompression
    }

    public enum ChemicalEquilibriumTransformType
    {
        TP,
        TV,
        HP,
        UV
    }

    public enum MixtureCompositionUpdateType
    {
        HP,
        UV
    }

    public enum CompressionType
    {
        Isentropic,
        Shock
    }

    public partial class IdealGasMixture
    {
        #region Private Fields
        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        [NonSerialized]
        IdealGasPotentialConservationLawSolver m_yuTemperatureSolver;

        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        [NonSerialized]
        IdealGasPotentialConservationLawSolver m_yhTemperatureSolver;

        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        [NonSerialized]
        IdealGasMolarEntropyConservationLawTemperatureSolver m_molarEntropyConservationLawTemperatureSolver;

        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        [NonSerialized]
        IdealGasMolarEntropyConservationLawPressureCalculator m_molarEntropyConservationLawPressureCalculator;

        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        [NonSerialized]
        IdealGasIsentropicCompressionTemperatureSolver m_isentropicCompressionTemperatureSolver;

        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        [NonSerialized]
        ShockWaveSoneSolver m_shockWaveSoneSolver; 
        #endregion

        StateChangeCategory EvaluateTPChangeCategory(double t, double p)
        {
            StateChangeCategory category = StateChangeCategory.None;
            if (false == m_temperature.UpToDate || m_temperature.Value != t)
            {
                category |= StateChangeCategory.Temperature;
            }
            if (false == m_pressure.UpToDate || m_pressure.Value != t)
            {
                category |= StateChangeCategory.Pressure;
            }
            return category;
        }

        public void Update(Temperature temperature, TemperatureChangeType changeType)
        {
            this.UpdateTemperature(temperature.ToDouble(), changeType);
        }

        [DebuggerStepThrough()]
        public void Update(Pressure pressure, PressureChangeType changeType)
        {
            this.UpdatePressure(pressure.ToDouble(), changeType);
        }

        public void Update(MassConcentration massDensity, DensityChangeType changeType)
        {
            this.UpdateMolarDensity(massDensity, changeType);
        }

        public void Update(CompressionRatio compressionRatio, CompressionType compressionType)
        {
            if (compressionRatio.Value == 1.0) { return; }
            switch (compressionType)
            {
                case (CompressionType.Isentropic):
                    {
                        var solver = m_isentropicCompressionTemperatureSolver ??
                            (m_isentropicCompressionTemperatureSolver = new IdealGasIsentropicCompressionTemperatureSolver(Model));
                        solver.Function.Configure(MolarS, this.Temperature, compressionRatio.Value, SpeciesMoleFractions);
                        double t = solver.Solve(this.Temperature);
                        double p = compressionRatio.Value * this.Pressure * t / this.Temperature;
                        UpdateTP(t, p);
                    }
                    break;
                case (CompressionType.Shock):
                    {
                        if (compressionRatio.Value < 1.0d)
                        {
                            throw new ArgumentOutOfRangeException("compressionRatio");
                        }
                        double gamma = this.Cp / this.Cv;
                        double crMax = (gamma + 1d) / (gamma - 1d);
                        DoubleRelativeTolerance tol = new DoubleRelativeTolerance(5.0e-2);
                        if( compressionRatio.Value > crMax ||
                            tol.AreSame(compressionRatio.Value, crMax))
                        {
                            double t = this.Temperature;
                            double p = this.Pressure;
                            double density = MassDensity;
                            UpdateTemperature(10000.0d, TemperatureChangeType.ShockCompression);
                            crMax = MassDensity / density;
                            if (compressionRatio.Value > crMax)
                            {
                                throw new ArgumentOutOfRangeException("compressionRatio",
                                        string.Format("The compression ratio must not be higher than {0}", crMax)
                                    );
                            }
                            UpdateTP(t, p);
                        }
                        var ratios = PolytropicGasRatios.CaclShockWaveRatios(new CompressionRatio(compressionRatio), gamma);
                        double[] x = new double[4];
                        var v = (ShockCompressionVector)x;
                        v.Temperature = this.Temperature * ratios.TemperatureRatio; 
                        v.Pressure = this.Pressure * ratios.PressureRatio;
                        v.CompressionRatio = ratios.DensityRatio;
                        v.MachNumber = Math.Sqrt(ratios.TemperatureRatio);

                        var solver = m_shockWaveSoneSolver ??
                            (m_shockWaveSoneSolver = new ShockWaveSoneSolver());
                        solver.Function.Configure(this, null, null, compressionRatio.Value, null);
                        object anchor = null;
                        solver.Solve(x, ref anchor);
                        UpdateTP(v.Temperature, v.Pressure);
                    }
                    break;
                default:
                    throw new NotSupportedException(compressionType.ToString());
            }
        }

        public void Update(MachNumber machNumber)
        {
            if (machNumber.Value < 1d)
            {
                throw new ArgumentOutOfRangeException("machNumber");
            }
            else if (machNumber.Value == 1d) { return; }
            double gamma = this.Cp / this.Cv;
            var ratios = PolytropicGasRatios.CaclShockWaveRatios(machNumber, gamma);
            double[] x = new double[4];
            var v = (ShockCompressionVector)x;
            v.Temperature = this.Temperature * ratios.TemperatureRatio;
            v.Pressure = this.Pressure * ratios.PressureRatio;
            v.CompressionRatio = ratios.DensityRatio;
            v.MachNumber = Math.Sqrt(ratios.TemperatureRatio);

            var solver = m_shockWaveSoneSolver ??
                (m_shockWaveSoneSolver = new ShockWaveSoneSolver());
            solver.Function.Configure(this, null, null, null, machNumber.Value);
            object anchor = null;
            solver.Solve(x, ref anchor);
            UpdateTP(v.Temperature, v.Pressure);
        }

        public void UpdateY(decimal[] y, MixtureCompositionUpdateType updateType)
        {
            using (var item = Model.SpeciesArraysPool.Acquire())
            {
                double[] yd = item.Unwrap();
                for (int isp = 0; isp < Nsp; ++isp)
                {
                    yd[isp] = (double)y[isp];
                }
                UpdateY(yd, updateType);
            }                 
        }

        public void UpdateX(decimal[] x, MixtureCompositionUpdateType updateType)
        {
            using (var item = Model.SpeciesArraysPool.Acquire())
            {
                double[] xd = item.Unwrap();
                for (int isp = 0; isp < Nsp; ++isp)
                {
                    xd[isp] = (double)x[isp];
                }
                UpdateX(xd, updateType);
            } 
        }

        public void UpdateY(double[] y, MixtureCompositionUpdateType updateType)
        {
            switch (updateType)
            {
                case (MixtureCompositionUpdateType.UV):
                    {
                        var solver = m_yuTemperatureSolver ??
                            (m_yuTemperatureSolver = IdealGasPotentialConservationLawSolver
                            .CreateSpecificInternalEnergyConservationLawTemperatureSolver(Model));
                        solver.Function.Configure(SpecificU, y);
                        double temperature = solver.Solve(this.Temperature);
                        UpdateTRO(temperature, this.MassDensity);
                    }
                    break;
                case (MixtureCompositionUpdateType.HP):
                    {
                        var solver = m_yhTemperatureSolver ??
                            (m_yhTemperatureSolver = IdealGasPotentialConservationLawSolver
                            .CreateSpecificEnthalpyConservationLawTemperatureSolver(Model));
                        solver.Function.Configure(SpecificH, y);
                        double temperature = solver.Solve(this.Temperature);
                        UpdateTP(temperature, this.Pressure);
                    }
                    break;
                default:
                    throw new NotSupportedException(updateType.ToString());
            }
        }

        public void UpdateX(double[] x, MixtureCompositionUpdateType updateType)
        {
            switch (updateType)
            {
                case (MixtureCompositionUpdateType.UV):
                    {
                        var solver = m_yuTemperatureSolver ??
                            (m_yuTemperatureSolver = IdealGasPotentialConservationLawSolver
                            .CreateMolarInternalEnergyConservationLawTemperatureSolver(Model));
                        solver.Function.Configure(MolarU, x);
                        double temperature = solver.Solve(this.Temperature);
                        UpdateTRO(temperature, this.MassDensity);
                    }
                    break;
                case (MixtureCompositionUpdateType.HP):
                    {
                        var solver = m_yhTemperatureSolver ??
                            (m_yhTemperatureSolver = IdealGasPotentialConservationLawSolver
                            .CreateMolarEnthalpyConservationLawTemperatureSolver(Model));
                        solver.Function.Configure(MolarH, x);
                        double temperature = solver.Solve(this.Temperature);
                        UpdateTP(temperature, this.Pressure);
                    }
                    break;
                default:
                    throw new NotSupportedException(updateType.ToString());
            }
        }

        public void UpdateTemperature(double t, TemperatureChangeType changeType)
        {
            if (t == this.Temperature) { return; }
            switch (changeType)
            {
                case( TemperatureChangeType.AtConstantPressure):
                    double p = m_pressure.Value;
                    DropState(StateChangeCategory.Temperature);
                    m_temperature.Value = t;
                    m_pressure.Value = p;
                    m_molarDensity.Value = p / t / Constants.Rgas;
                    break;
                case (TemperatureChangeType.AtConstantVolume):
                    double c = m_molarDensity.Value;
                    DropState(StateChangeCategory.Temperature);
                    m_temperature.Value = t;                    
                    m_molarDensity.Value = c;
                    m_pressure.Value = c* Constants.Rgas * t;
                    break;
                case(TemperatureChangeType.IsentropicCompression):
                    {
                        var calculator = m_molarEntropyConservationLawPressureCalculator ??
                            (m_molarEntropyConservationLawPressureCalculator = new IdealGasMolarEntropyConservationLawPressureCalculator(Model));
                        double finalPressure = calculator.CalcPressure(t, this.Pressure, this.MolarS, this.SpeciesMoleFractions, this.SpeciesMoleFractionsLog);
                        UpdateTP(t, finalPressure);
                    }
                    break;
                case (TemperatureChangeType.ShockCompression):
                    {
                        var ratios = PolytropicGasRatios.CaclShockWaveRatios(new Temperature(this.Temperature), new Temperature(t), this.Cp / this.Cv);
                        double[] x = new double[4];
                        var v = (ShockCompressionVector)x;
                        v.Temperature = this.Temperature * ratios.TemperatureRatio; 
                        v.Pressure = this.Pressure * ratios.PressureRatio;
                        v.CompressionRatio = ratios.DensityRatio;
                        v.MachNumber = Math.Sqrt(ratios.TemperatureRatio);
                        var solver = m_shockWaveSoneSolver ??
                            (m_shockWaveSoneSolver = new ShockWaveSoneSolver());
                        solver.Function.Configure(this, t, null, null, null);
                        object anchor = new object();
                        solver.Solve(x, ref anchor);
                        this.UpdateTP(v.Temperature, v.Pressure);
                    }
                    break;
                default:
                    throw new NotSupportedException(changeType.ToString());
            }            
        }

        public void UpdatePressure(double p, PressureChangeType changeType)
        {
            if (p == this.Pressure) { return; }
            switch (changeType)
            {
                case( PressureChangeType.AtConstantTemperature ):
                    double t = m_temperature.Value;
                    DropState(StateChangeCategory.Pressure | StateChangeCategory.Density);
                    m_temperature.Value = t;
                    m_pressure.Value = p;
                    m_molarDensity.Value = p / t / Constants.Rgas;
                    break;
                case( PressureChangeType.AtConstantVolume):
                    double c = m_molarDensity.Value;
                    DropState(StateChangeCategory.Temperature | StateChangeCategory.Pressure);
                    m_molarDensity.Value = c;
                    m_pressure.Value = p;
                    m_temperature.Value = p / c / Constants.Rgas;
                    break;
                case (PressureChangeType.IsentropicCompression):                    
                    {
                        var solver = m_molarEntropyConservationLawTemperatureSolver ??
                            (m_molarEntropyConservationLawTemperatureSolver = new IdealGasMolarEntropyConservationLawTemperatureSolver(Model));
                        solver.Function.Configure(MolarS, p, SpeciesMoleFractions);
                        UpdateTP(solver.Solve(Temperature), p);
                    }
                    break;
                case (PressureChangeType.ShockCompression):
                    {
                        var ratios = PolytropicGasRatios.CaclShockWaveRatios(new Pressure(this.Pressure), new Pressure(p), this.Cp / this.Cv);
                        double[] x = new double[4];
                        var v = (ShockCompressionVector)x;
                        v.Temperature = this.Temperature * ratios.TemperatureRatio; ;
                        v.Pressure = this.Pressure * ratios.PressureRatio; ;
                        v.CompressionRatio = ratios.DensityRatio;
                        v.MachNumber = Math.Sqrt(ratios.TemperatureRatio);
                        var solver = m_shockWaveSoneSolver ??
                            (m_shockWaveSoneSolver = new ShockWaveSoneSolver());
                        solver.Function.Configure(this, null, p, null, null);                                                
                        object anchor = new object();
                        solver.Solve(x, ref anchor);
                        this.UpdateTP(v.Temperature, v.Pressure);
                    }
                    break;
                default:
                    throw new NotSupportedException(changeType.ToString());
            }  
        }

        public void UpdateMolarDensity(double molarDensity, DensityChangeType changeType)
        {
            switch (changeType)
            {
                case(DensityChangeType.AtConstantPressure):
                    double t = this.Pressure / molarDensity / Constants.Rgas;
                    UpdateTP(t, this.Pressure);
                    break;
                case (DensityChangeType.AtConstantTemperature):
                    double p = this.Temperature * Constants.Rgas * molarDensity;
                    UpdateTP(p, this.Temperature);
                    break;
                case (DensityChangeType.IsentropicCompression):
                    Update(new CompressionRatio(molarDensity / this.MolarDensity), CompressionType.Isentropic);
                    break;
                case (DensityChangeType.ShockCompression):
                    Update(new CompressionRatio(molarDensity / this.MolarDensity), CompressionType.Shock);
                    break;
                default:
                    throw new NotSupportedException(changeType.ToString());
            }
        }

        public void Update(IMixtureDefinition mixtureDefinition)
        {
            DropState(StateChangeCategory.Temperature | StateChangeCategory.Pressure | StateChangeCategory.Composition);
            m_temperature.Value = mixtureDefinition.Temperature;
            m_pressure.Value = mixtureDefinition.Pressure;
            m_molarDensity.Value = m_pressure.Value / m_temperature.Value / Constants.Rgas;
            double[] x = new double[Nsp];
            for (int jsp = 0; jsp < mixtureDefinition.X.Length; ++jsp)
            {
                int isp = Model.GetSpeciesIndex( mixtureDefinition.SpeciesIDs[jsp]);
                x[isp] = mixtureDefinition.X[jsp];
            }
            m_speciesMoleFractions.Assign(x);
        }

        public void Update(ChemicalEquilibriumTransformType changeType)
        {
            using (var xitem = Model.SpeciesArraysPool.Acquire())
            {
                double[] x = xitem.Unwrap();
                this.SpeciesMoleFractions.CopyTo(x);
                double t = this.Temperature;
                double p = this.Pressure;
                switch (changeType) 
                {
                    case(ChemicalEquilibriumTransformType.TP):
                        Model.IdealGasChemicalEquilibriumCalculator.TP(t, p, x);
                        UpdateTPX(t, p, x);
                        break;
                    case (ChemicalEquilibriumTransformType.TV):   
                        Model.IdealGasChemicalEquilibriumCalculator.TV(t, ref p, x);
                        UpdateTPX(t, p, x);
                        break;
                    case (ChemicalEquilibriumTransformType.HP):
                        Model.IdealGasChemicalEquilibriumCalculator.HP(ref t, p, x);
                        UpdateTPX(t, p, x);
                        break;
                    case (ChemicalEquilibriumTransformType.UV):
                        Model.IdealGasChemicalEquilibriumCalculator.UV(ref t, ref p, x);
                        UpdateTPX(t, p, x);
                        break;
                    default:
                        throw new NotSupportedException(changeType.ToString());
                }
            }
        }

        public void UpdateTP(double t, double p)
        {
            var category = EvaluateTPChangeCategory(t, p);
            DropState(category);
            m_temperature.Value = t;
            m_pressure.Value = p;
            m_molarDensity.Value = m_pressure.Value / m_temperature.Value / Constants.Rgas;
        }

        public void UpdateTRO(double t, double ro)
        {
            double p = Constants.Rgas * t * ro / this.MolarMass;
            UpdateTP(t, p);
        }

        public void UpdateTPX(double t, double p, ReadOnlyArray<double> x)
        {
            StateChangeCategory category = EvaluateTPChangeCategory(t, p) | StateChangeCategory.Composition;
            DropState(category);
            m_temperature.Value = t;
            m_pressure.Value = p;
            m_molarDensity.Value = p / t / Constants.Rgas;
            m_speciesMoleFractions.Assign(x.Data);
        }

        public void UpdateTPY(double t, double p, ReadOnlyArray<double> y)
        {
            StateChangeCategory category = EvaluateTPChangeCategory(t, p) | StateChangeCategory.Composition;
            DropState(category);
            m_temperature.Value = t;
            m_pressure.Value = p;
            m_molarDensity.Value = p / t / Constants.Rgas;
            m_speciesMassFractions.Assign(y.Data);
        }

        public void UpdateTPX(IdealGasMixture other)
        {
            this.UpdateTPX(other.Temperature, other.Pressure, other.SpeciesMoleFractions);
        }

        public void UpdateTPY(IdealGasMixture other)
        {
            this.UpdateTPX(other.Temperature, other.Pressure, other.SpeciesMassFractions);
        }

        public void IncrementX(ReadOnlyArray<double> dx, MixtureCompositionUpdateType updateType)
        {
            using( var item = Model.SpeciesArraysPool.Acquire())
            {
                double[] xnew = item.Unwrap();
                double[] xold = this.SpeciesMoleFractions.Data;
                for (int isp = 0; isp < Nsp; ++isp)
                {
                    xnew[isp] = xold[isp] + dx[isp];
                }
                UpdateX(xnew, updateType);
            }
        }

        public void IncrementY(ReadOnlyArray<double> dy, MixtureCompositionUpdateType updateType)
        {
            using (var item = Model.SpeciesArraysPool.Acquire())
            {
                double[] ynew = item.Unwrap();
                double[] yold = this.SpeciesMoleFractions.Data;
                for (int isp = 0; isp < Nsp; ++isp)
                {
                    ynew[isp] = yold[isp] + dy[isp];
                }
                UpdateY(ynew, updateType);
            }
        }

        public void IncrementX(ReadOnlyArray<decimal> dx, MixtureCompositionUpdateType updateType)
        {
            using (var item = Model.SpeciesArraysPool.Acquire())
            {
                double[] xnew = item.Unwrap();
                double[] xold = this.SpeciesMoleFractions.Data;
                for (int isp = 0; isp < Nsp; ++isp)
                {
                    xnew[isp] = xold[isp] + (double)dx[isp];
                }
                UpdateX(xnew, updateType);
            }
        }

        public void IncrementY(ReadOnlyArray<decimal> dy, MixtureCompositionUpdateType updateType)
        {
            using (var item = Model.SpeciesArraysPool.Acquire())
            {
                double[] ynew = item.Unwrap();
                double[] yold = this.SpeciesMoleFractions.Data;
                for (int isp = 0; isp < Nsp; ++isp)
                {
                    ynew[isp] = yold[isp] + (double)dy[isp];
                }
                UpdateY(ynew, updateType);
            }
        }

    }

    class ShockWaveSoneFunction : ISoneFunction
    {
        #region Private Fields
        IdealGasMixture m_idealGasMixture;
        double m_h0, m_t0, m_ro0, m_p0, m_c0;
        double m_fixedCompressionRatio;
        double m_fixedMachNumber;
        double m_fixedTemperature;
        double m_fixedPressure;
        double m_fixedEnthalpy;
        Action<double[], double[]> m_selectedFormulation;
        #endregion

        internal void Configure(IdealGasMixture idealGasMixture, 
            double? temperature, 
            double? pressure,
            double? compressionRatio,
            double? machNumber)
        {
            m_idealGasMixture = idealGasMixture;
            m_t0 = idealGasMixture.Temperature;
            m_p0 = idealGasMixture.Pressure;
            m_ro0 = idealGasMixture.MassDensity;
            m_h0 = idealGasMixture.MolarH;
            m_c0 = idealGasMixture.SpeedOfSound;

            if (compressionRatio.HasValue)
            {
                m_fixedCompressionRatio = compressionRatio.Value;
                m_selectedFormulation = this.CompressionRatioFormulation;
            }
            else if (machNumber.HasValue)
            {
                m_fixedMachNumber = machNumber.Value;
                m_selectedFormulation = this.MachNumberFormulation;
            }
            else if (pressure.HasValue && !temperature.HasValue)
            {
                m_fixedPressure = pressure.Value;
                m_selectedFormulation = this.FixedPressureFormulation;
            }
            else if (temperature.HasValue && !pressure.HasValue)
            {
                m_idealGasMixture.UpdateTemperature(temperature.Value, TemperatureChangeType.AtConstantVolume);
                m_fixedTemperature = temperature.Value;
                m_fixedEnthalpy = idealGasMixture.MolarH;
                m_selectedFormulation = this.FixedTemperatureFormulation;
            }
            else
            {
                throw new ArgumentException();
            }
        }

        #region ISoneFunction Members

        public int Nva
        {
            get { return 4; }
        }

        public void Compute(double[] x, double[] y)
        {
            m_selectedFormulation(x, y);
        }

        private void MachNumberFormulation(double[] x, double[] y)
        {
            ShockCompressionVector xv = (ShockCompressionVector)x;
            m_idealGasMixture.UpdateTP(xv.Temperature, xv.Pressure);

            double t = xv.Temperature;
            double p = xv.Pressure;
            double m = xv.MachNumber;
            double cr = xv.CompressionRatio;
            double h = m_idealGasMixture.MolarH;
            double ro = m_idealGasMixture.MassDensity;
            double d = xv.MachNumber * m_c0;
            double j = m * m_ro0 * m_c0;

            y[0] = (m_h0 - h) + 0.5 * Constants.Rgas * (t / p + m_t0 / m_p0) * (p - m_p0);
            y[1] = j * j - (p - m_p0) / (1.0 / m_ro0 - 1.0 / ro);
            y[2] = cr - ro / m_ro0;
            y[3] = m - m_fixedMachNumber;
        }


        private void CompressionRatioFormulation(double[] x, double[] y)
        {
            ShockCompressionVector xv = (ShockCompressionVector)x;
            m_idealGasMixture.UpdateTP(xv.Temperature, xv.Pressure);

            double t = xv.Temperature;
            double p = xv.Pressure;
            double m = xv.MachNumber;
            double cr = xv.CompressionRatio;
            double h = m_idealGasMixture.MolarH;
            double d = xv.MachNumber * m_c0;
            double j = m * m_ro0 * m_c0;
            double w = m_idealGasMixture.MolarMass;

            double fcr = m_fixedCompressionRatio;

            y[0] = (m_h0 - h) / w + 0.5 * m_p0 / m_ro0 * (fcr * t / m_t0 - 1.0) * (fcr + 1.0) / fcr;
            y[1] = m - 1.0 / m_ro0 / m_c0 *
                Math.Sqrt(
                    fcr * m_p0 * m_ro0 * (fcr * t / m_t0 - 1.0) / (fcr - 1.0)
                );
            y[2] = p - m_p0 * fcr * t / m_t0;
            y[3] = cr - fcr;
        }

        private void FixedTemperatureFormulation(double[] x, double[] y)
        {
            ShockCompressionVector xv = (ShockCompressionVector)x;
            m_idealGasMixture.UpdateTP(xv.Temperature, xv.Pressure);

            double fxt = m_fixedTemperature;
            double t = xv.Temperature;
            double p = xv.Pressure;
            double m = xv.MachNumber;
            double cr = xv.CompressionRatio;
            double h = m_fixedEnthalpy;
            double ro = m_idealGasMixture.MassDensity;
            double d = xv.MachNumber * m_c0;
            double j = m * m_ro0 * m_c0;
            double w = m_idealGasMixture.MolarMass;

            y[0] = (m_h0 - h) / w + 0.5 * m_p0 / m_ro0 * (cr * fxt / m_t0 - 1.0) * (cr + 1.0) / cr;
            y[1] = m - 1.0 / m_ro0 / m_c0 *
                Math.Sqrt(
                    cr * m_p0 * m_ro0 * (cr * fxt / m_t0 - 1.0) / (cr - 1.0)
                );
            y[2] = p - m_p0 * cr * fxt / m_t0;
            y[3] = t - fxt;
        }


        private void FixedPressureFormulation(double[] x, double[] y)
        {
            ShockCompressionVector xv = (ShockCompressionVector)x;
            m_idealGasMixture.UpdateTP(xv.Temperature, xv.Pressure);
            double fxp = m_fixedPressure;

            double t = xv.Temperature;
            double p = xv.Pressure;
            double m = xv.MachNumber;
            double cr = xv.CompressionRatio;
            double h = m_idealGasMixture.MolarH;
            double ro = m_idealGasMixture.MassDensity;
            double d = xv.MachNumber * m_c0;
            double j = m * m_ro0 * m_c0;
            double w = m_idealGasMixture.MolarMass;


            y[0] = (m_h0 - h) + 0.5 * (fxp - m_p0) * Constants.Rgas * (t / fxp + m_t0 / m_p0);
            y[1] = m - 1.0 / m_ro0 / m_c0 *
                Math.Sqrt(
                        (fxp - m_p0) / (1.0 / m_ro0 - 1.0 / ro)
                    );
            y[2] = cr - ro / m_ro0;
            y[3] = p - fxp;
        }

        #endregion
    }

    class ShockWaveSoneSolver : SoneNewton<ShockWaveSoneFunction>
    {
        internal ShockWaveSoneSolver()
            : base(new ShockWaveSoneFunction())
        {
            SetAbsoluteConvergenceRadius(
                    ShockCompressionVector.Create(100.0, 10.0 * Constants.Atmosphere, 10.0, 0.5)
                );
            SetRelativeConvergenceRadius(
                    ShockCompressionVector.Create(1.0e-1, 0.5, 5.0, 0.5)
                );

            SetAbsoluteTolerance(0.0);
            SetRelativeTolerance(1.0e-8);
        }



    }

    /// <summary>
    /// ShockCompressionVector
    /// </summary>
    partial struct ShockCompressionVector : Nature.Numerics.IVector
    {
        #region Indexes


        /// <summary>
        /// Temperature behind the shock wave
        /// </summary>
        public const int TemperatureIndex = 0;


        /// <summary>
        /// Pressure behind the shock wave
        /// </summary>
        public const int PressureIndex = 1;


        /// <summary>
        /// Compression Ratio
        /// </summary>
        public const int CompressionRatioIndex = 2;


        /// <summary>
        /// Mach number: U / C
        /// </summary>
        public const int MachNumberIndex = 3;


        #endregion

        #region Properties


        /// <summary>
        /// Temperature behind the shock wave
        /// </summary>
        public double Temperature
        {
            get { return m_array[0]; }
            set { m_array[0] = value; }
        }


        /// <summary>
        /// Pressure behind the shock wave
        /// </summary>
        public double Pressure
        {
            get { return m_array[1]; }
            set { m_array[1] = value; }
        }


        /// <summary>
        /// Compression Ratio
        /// </summary>
        public double CompressionRatio
        {
            get { return m_array[2]; }
            set { m_array[2] = value; }
        }


        /// <summary>
        /// Mach number: U / C
        /// </summary>
        public double MachNumber
        {
            get { return m_array[3]; }
            set { m_array[3] = value; }
        }


        #endregion

        #region Constructors

        public static ShockCompressionVector Create()
        {
            ShockCompressionVector vector = new ShockCompressionVector();
            vector.m_array = new double[4];
            return vector;
        }

        public static ShockCompressionVector Create(double temperature, double pressure, double compressionRatio, double machNumber)
        {
            ShockCompressionVector vector = new ShockCompressionVector();
            vector.m_array = new double[4];

            vector.m_array[0] = temperature;


            vector.m_array[1] = pressure;


            vector.m_array[2] = compressionRatio;


            vector.m_array[3] = machNumber;


            return vector;
        }

        private ShockCompressionVector(double[] array)
            : this()
        {
            if (array.Length != 4) { throw new ArgumentOutOfRangeException(); }
            m_array = array;
        }

        public static explicit operator ShockCompressionVector(double[] array)
        {
            return new ShockCompressionVector(array);
        }

        public static implicit operator double[](ShockCompressionVector vector)
        {
            return vector.m_array;
        }

        #endregion

        #region Private Fields

        [DebuggerBrowsable(DebuggerBrowsableState.RootHidden)]
        private double[] m_array;

        #endregion

        #region IVector Members

        public double this[int i]
        {
            get { return m_array[i]; }
            set { m_array[i] = value; }
        }

        public int Length
        {
            get { return m_array.Length; }
        }

        #endregion

        #region IEnumerable<double> Members

        public System.Collections.Generic.IEnumerator<double> GetEnumerator()
        {
            for (int i = 0; i < m_array.Length; ++i) { yield return m_array[i]; }
        }

        #endregion

        #region IEnumerable Members

        System.Collections.IEnumerator System.Collections.IEnumerable.GetEnumerator()
        {
            return m_array.GetEnumerator();
        }

        #endregion

    }
}
